!>\brief
!! In this module the local matrix and right hand side for the various
!! semi-conductor problems are computed.
!!
!! \n
!!
!! The local matrix does not include any static condensation; it
!! includes the Neumann type boundary conditions but not the Dirichlet
!! type ones.
!<----------------------------------------------------------------------
module mod_sc_locmat

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me

 use mod_grid, only: &
   mod_grid_initialized, &
   t_e, affmap

 use mod_bcs, only: &
   mod_bcs_initialized, &
   p_t_b_v, p_t_b_s, p_t_b_e, &
   b_dir,   b_neu,   b_ddc

 use mod_scharfetter_gummel, only: &
   mod_scharfetter_gummel_initialized, &
   bim3a_local_laplacian, &
   bim3a_osc_local_laplacian, &
   bim3a_local_scharfetter_gummel, &
   bim3a_osc_local_scharfetter_gummel, &
   bim3a_local_reaction, &
   bim3a_local_rhs

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_sc_locmat_constructor, &
   mod_sc_locmat_destructor,  &
   mod_sc_locmat_initialized, &
   sc_locmat, &
   sc_masmat, &
   t_bcs_error

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! This type is used to indicate that an error has occurred in the
 ! computation of the boundary condition
 type t_bcs_error
   logical :: lerr
   character(len=1000) :: message
 end type t_bcs_error

! Module variables

 ! public members
 logical, protected ::               &
   mod_sc_locmat_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_sc_locmat'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_sc_locmat_constructor()

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_bcs_initialized.eqv..false.)      .or. &
(mod_scharfetter_gummel_initialized.eqv..false.))then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_sc_locmat_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_sc_locmat_initialized = .true.
 end subroutine mod_sc_locmat_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_sc_locmat_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_sc_locmat_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_sc_locmat_initialized = .false.
 end subroutine mod_sc_locmat_destructor

!-----------------------------------------------------------------------

 !> Computation of the local mass matrix
 pure subroutine sc_masmat( mmk , e )
  !> local matrix and rhs
  real(wp),      intent(out) :: mmk(:,:)
  type(t_e),     intent(in) :: e  !< element

  real(wp) :: node_coeff(size(mmk,1))
 
   node_coeff = 1.0_wp
   call bim3a_local_reaction( e%vol,node_coeff,1.0_wp,mmk )

 end subroutine sc_masmat

!-----------------------------------------------------------------------
 
 !> Computation of the local matrix and right hand side.
 !!
 !! Since we assume \f$\mathbb{P}_1\f$ finite elements, the gradients
 !! in physical space are readily obtained by
 !! \f{displaymath}{
 !!  \nabla\varphi_i = - \frac{1}{h_i} \underline{n}_i
 !! \f}
 !! where \f$\underline{n}_i\f$ is the outward normal of the
 !! \f$i\f$-side (i.e. the side
 !! opposed to the \f$i\f$-vertex, see \c mod_grid) and \f$h_i\f$ is
 !! the \f$i\f$-th height, obtained as
 !! \f{displaymath}{
 !!  h_i = d\frac{|K|}{|e_i|}.
 !! \f}
 !!
 !! The input values refer to the following problem:
 !! \f{displaymath}{
 !!  -\nabla\cdot\left( \tilde{\mu}\left( \nabla u +
 !!  u\nabla\tilde{\varphi} \right) \right) = f.
 !! \f}
 !! This implies that for electron and holes densities we need to pass
 !! as \c mu_e and \c v_v the scaled mobility and potential
 !! \f$\mu_{p,n}V_{th}\f$ and \f$\frac{\varphi}{V_{th}}\f$,
 !! respectively.
 !!
 !! \note If the input argument \c v_v is not present, it is assumed
 !! that there is no transport term.
 !! \note Zero-order terms are not considered since their matrix is
 !! diagonal and they are better treated at the global level.
 !!
 !! The detailed documentation concerning the computation of the local
 !! matrices is given in \c mod_scharfetter_gummel, in the following
 !! procedures:
 !! <ul>
 !!  <li> \c bim3a_local_scharfetter_gummel: stabilized advection
 !!  diffusion matrix
 !!  <li> \c bim3a_local_reaction: reaction matrix
 !!  <li> \c bim3a_local_rhs: right hand side.
 !! </ul>
 !!
 !! \warning In this subroutine, the nodal and element right-hand side
 !! contributions \c f_v and \c f_e are <em>added</em>, which is
 !! different from waht it done in \c bim3a_local_rhs where such
 !! contributions are intended to be <em>multiplied</em>. In practice,
 !! this means that here we only pass nodal values to bim3a_local_rhs.
 pure subroutine sc_locmat( aak, fk, me, e, be, bcs_err, &
                            mu_e, f_v, f_e, v_v )
  !> local matrix and rhs
  real(wp),      intent(out) :: aak(:,:), fk(:)
  type(t_me),    intent(in) :: me !< master element
  type(t_e),     intent(in) :: e  !< element
  type(p_t_b_e), intent(in) :: be !< boundary element
  real(wp),      intent(in) :: mu_e !> element diffusion coeff.
  !> right hand side: nodal values
  real(wp),      intent(in), optional :: f_v(:)
  !> right hand side: element value
  real(wp),      intent(in), optional :: f_e
  !> (scaled) nodal potential (present only for advection-diffusion
  !! problems)
  real(wp),      intent(in), optional :: v_v(:)
  type(t_bcs_error), intent(out) :: bcs_err !< bcs error structure
 
  integer :: iv
  real(wp) :: gradp(e%d,e%d+1), xv(e%d,e%d+1), f_ve(e%d+1)
  character(len=*), parameter :: &
    this_sub_name = 'sc_locmat'

   ! compute the gradients of the basis functions in physical space
   do iv=1,e%d+1
      xv(:,iv) = e%v(iv)%p%x
     gradp(:,iv) = -e%s(iv)%p%a/(real(e%d,wp)*e%vol) * e%n(:,iv)
   enddo

   if(present(v_v)) then
      !call bim3a_local_scharfetter_gummel (me, gradp, e%vol, mu_e, v_v, aak)
      call bim3a_osc_local_scharfetter_gummel (me, xv, gradp, e%vol, mu_e, v_v, aak)
   else
      !call bim3a_local_laplacian (gradp, e%vol, mu_e, aak)
      call bim3a_osc_local_laplacian (xv, gradp, e%vol, mu_e, aak)
   endif

   f_ve = 0.0_wp
   if(present(f_v)) f_ve = f_ve + f_v ! nodal values
   if(present(f_e)) f_ve = f_ve + f_e ! constant value for all the nodes
   call bim3a_local_rhs (e%vol, f_ve, 1.0_wp, fk)

   ! add boundary fluxes if present
   bcs_err%lerr = .false.
!   if(associated(be%p)) then
!     call sc_boundary_flux(aak,fk,base,be,bcs_err)
!   endif
   !--------------------------------------------------------------------

 end subroutine sc_locmat
 
!-----------------------------------------------------------------------
 
! pure subroutine sc_boundary_flux(aak,fk,base,be,err)
!  real(wp), intent(inout) :: aak(:,:), fk(:)
!  type(t_base), intent(in) :: base
!  type(p_t_b_e), intent(in) :: be
!  type(t_bcs_error), intent(out) :: err
!
!  integer :: isb, isl, breg, l, i, j
!                              !integer l_e(base%me)
!  real(wp) :: xgs(be%p%e%m,base%ms), a(be%p%e%m,base%ms), &
!    gamma(base%ms), hb(base%ms), an(base%ms), ap(base%ms),&
!    wgsa(base%ms)
! 
!   err%lerr = .false.
!   do isb=1,be%p%nbs ! loop on the boundary sides
!     btype: select case(be%p%bs(isb)%p%bc)
!      case(b_dir,b_ddc)
!       ! nothing to do in this subroutine
!      case(b_neu)
! 
!       isl = be%p%bs(isb)%p%s%isl(1) ! local index of the side
!       ! Gaussian nodes in physical space
!       xgs = affmap(be%p%e,base%xigb(:,:,isl))
!       ! side weigths reordered and scaled by the side area
!       wgsa = (be%p%e%s(isl)%p%a/base%me%voldm1) * &
!                   base%wgs(base%stab(be%p%e%ip(isl),:))
!
!       ! problem coefficients
!       breg = -be%p%bs(isb)%p%s%ie(2) ! boundary region
!       gamma = coeff_rob(xgs,breg)
!       hb    = coeff_neu(xgs,breg)
!
!       btypespec: select case(be%p%bs(isb)%p%btype)
!
!        case(1) ! Neumann, Hughes flux
!
!         ! advection coefficient
!         a  = coeff_adv(xgs)
!         an = matmul(transpose(a),be%p%e%n(:,isl))
!         ap = 0.5_wp*(an + abs(an)) + gamma
!
!         do l=1,base%ms
!           do i=1,base%pk
!             do j=1,base%pk
!               aak(i,j) = aak(i,j) + wgsa(l) * &
!                 base%pb(j,l,isl) * ap(l) * base%pb(i,l,isl)
!             enddo
!             fk(i) = fk(i) - wgsa(l) * hb(l) * base%pb(i,l,isl)
!           enddo
!         enddo
!         
!        case default
!         err%lerr = .true.
!         write(err%message,'(A,I5,A,I5,A,I5,A)') &
!           'Unknown boundary condition type "',be%p%bs(isb)%p%btype, &
!           '" on side ',be%p%bs(isb)%p%s%i,', on element ',be%p%e%i,'.'
!       end select btypespec
!      case default
!       err%lerr = .true.
!       write(err%message,'(A,I5,A,I5,A,I5,A)') &
!         'Unknown boundary condition "',be%p%bs(isb)%p%bc,         &
!         '" on side ',be%p%bs(isb)%p%s%i,', on element ',be%p%e%i,'.'
!     end select btype
!   enddo

! end subroutine sc_boundary_flux
 
!-----------------------------------------------------------------------

end module mod_sc_locmat

